Subroutine calsimeq(xx,pres,phi,totmol)

	USE ckvar

	REAL(dp) :: fj(nspc),matj(nspc,nspc),df(nspc)
	REAL(dp) :: xx(nspc),xold(nspc)
	REAL :: pres
	REAL :: phi
	REAL :: errmx,err(nspc)
	REAL :: totmol,xtot
	REAL(dp) :: rxx
	REAL :: aa
	INTEGER :: i,j,k,n
	INTEGER :: nunk	! # of unknowns

	open(unit=60,file='out.dat')

	nunk=nspc
	aa=0.5/phi

	do 10 n=1,100
	xold=xx
	fj=0.0
	matj=0.0

! set up Kp equations for each reaction
! set up 1st derivation matrix of Kp equations
	do i=1,nreac
!	write(*,*) 'Reaction#:',i,'Kp=',kp(i)
	fj(i)=1.0
	rxx=kp(i)
	  do j=1,4
	  if(nspreac(i,j).gt.0.0) then
	    fj(i)=fj(i)*(xx(ispreac(i,j))*pres)**(nspreac(i,j))
	  elseif(nspreac(i,j).lt.0.0) then
	    rxx=rxx/(xx(ispreac(i,j))*pres)**(nspreac(i,j))
	  endif
	  enddo
	  do j=1,nunk
	    do k=1,4
	    if(j.eq.ispreac(i,k)) then
	      if(nspreac(i,k).gt.0.0) then
	        matj(i,j)=fj(i)*nspreac(i,k)/xx(j)
	      elseif(nspreac(i,k).lt.0.0) then
		matj(i,j)=rxx*nspreac(i,k)/xx(j)
	      endif
!	    write(*,*) spc(ispreac(i,k)),xx(j),matj(i,j)
	    endif
	    enddo
!	  write(*,*) j,matj(i,j)
	  enddo
	fj(i)=fj(i)-rxx
	enddo

!! set up atom conservations
! 1. Sum of every mole fraction:Eq.6
	do i=1,nunk
	fj(nreac+1)=fj(nreac+1)+xx(i)
	matj(nreac+1,i)=1.0
	enddo
	fj(nreac+1)=fj(nreac+1)-1.0
! 2. # of H / # of O: Eq.7
	fj(nreac+2)=2.0*aa*xx(1)-2.0*xx(2)+(2.0*aa-1.0)*xx(4)+&
		    (aa-1.0)*xx(5)-xx(6)+aa*xx(7)-xx(8)
	matj(nreac+2,1)=2.0*aa
	matj(nreac+2,2)=-2.0
	matj(nreac+2,3)=0.0
	matj(nreac+2,4)=2.0*aa-1.0
	matj(nreac+2,5)=aa-1.0
	matj(nreac+2,6)=-1.0
	matj(nreac+2,7)=aa
	matj(nreac+2,8)=-1.0
! 4. # of O / # of N: Eq.8 
	fj(nreac+3)=7.52*xx(2)-2.0*xx(3)+3.76*(xx(4)+xx(5)+xx(6))+&
		    2.76*xx(8)
	matj(nreac+3,1)=0.0
	matj(nreac+3,2)=7.52
	matj(nreac+3,3)=-2.0
	matj(nreac+3,4)=3.76
	matj(nreac+3,5)=3.76
	matj(nreac+3,6)=3.76
	matj(nreac+3,7)=0.0
	matj(nreac+3,8)=2.76
	
!	write(60,*) 'Iteration#:',n
	write(60,*) 'xx:'
	write(60,200) (xx(i),i=1,nspc)
	write(60,*) 'fj:'
	write(60,200) (fj(i),i=1,nspc)
	write(60,*) 'matj:'
	do i=1,nspc
	write(60,200) (matj(i,j),j=1,nspc)
	enddo

	CALL GAUSS(matj,nspc,df,-fj)

	write(60,*) 'df:'
	write(60,200) (df(j),j=1,nspc)
	write(60,*) 'xx+df'
	write(60,200) (xx(i)+df(i),i=1,nspc)

!	CALL LU(matj,df,-fj)

	errmx=0.0
	do i=1,nunk
	xx(i)=max(1.0e-30,xx(i)+df(i))
	err(i)=abs((xx(i)-xold(i))/xold(i))*100.0
	errmx=max(err(i),errmx)
	enddo
	write(60,*) 'Iteration#:',n,'errmx=',errmx
	if(errmx.le.0.001) go to 20
  10	continue
  20	continue
!	xtot=0.0
!	do i=1,nspc
!	xtot=xtot+xx(i)
!	enddo
!	write(*,*) 'errmx=',errmx,'xtot=',xtot
	totmol=7.52*aa/(2.0*xx(3)+xx(8))
  200	format(8(1pe9.2,1x))
End Subroutine

Subroutine LU(a,x,b)
! LU DECOMPOSITION MADE BY J.C. KIM :
	INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,100)

	REAL(dp) :: a(8,8),x(8),b(8)
	REAL(dp) :: factor,sum
	integer :: n,i,j,k
	n=8

	do k=1,n-1
	  do i=k+1,n
	    factor=a(i,k)/a(k,k)
	    a(i,k)=factor
	    do j=k+1,n
	      a(i,j)=a(i,j) - factor*a(k,j)
	    enddo
	  enddo
	enddo

	do i=2,n
	  sum=b(i)
	  do j=1,i-1
	    sum=sum - a(i,j)*b(j)
	  enddo
	  b(i)=sum
	enddo

	x(n)=b(n)/a(n,n)
	do i=n-1,1,-1
	  sum=0.0
	  do j=i+1,n
	    sum=sum+a(i,j)*x(j)
	  enddo
	  x(i)=(b(i)-sum)/a(i,i)
	enddo

	return
end subroutine

Subroutine GAUSS(a,n,dx,f)
	
	IMPLICIT NONE	
	INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,100)

	INTEGER :: n
	REAL(dp) :: a(n,n),dx(n),f(n)
	REAL(dp) :: b(n,n),c,d,tmp(n),sum
	INTEGER :: i,j,k,m,imax(1),ipvt(n)

	b=a
	ipvt=(/(i,i=1,n)/)
	
	do k=1,n
	imax=MAXLOC(ABS(b(k:n,k)))
	m=k-1+imax(1)

	if(m.ne.k) then
	ipvt((/m,k/))=ipvt((/k,m/))
	b((/m,k/),:)=b((/k,m/),:)
	endif
	d=1.0/b(k,k)

	tmp=b(:,k)
	  do j=1,n
	  c=b(k,j)*d
	  b(:,j)=b(:,j)-tmp*c
	  b(k,j)=c
	  enddo
	b(:,k)=tmp*(-d)
	b(k,k)=d
	enddo

	a(:,ipvt)=b

	do k=1,n
	sum=0.0
	  do j=1,n
	  sum=sum+a(k,j)*f(j)
	  enddo
	dx(k)=sum
	enddo
End subroutine
